/* 
 * bessel.cpp
 * 
 * Kirankumar R. Hiremath
 * Applied Analysis and Mathematical Physics Group
 * University of Twente, Department of Applied Mathematics 
 * P.O. Box 217, 7500AE Enschede, The Netherlands
 * (2004)
 *
 * email: k.r.hiremath@ieee.org
 */

/*
  Description: 
  1.This is a C++ library (part of it based on FORTRAN) to compute (double 
    precision) Bessel functions of first kind (J), bessel function of second 
    kind (Y), Hankel function of first kind (H1) and Hankel function of second 
    kind (H2); and their first derivative w.r.t. argument for REAL/COMPLEX 
    ORDER and REAL/COMPLEX ARGUMENT. 
  2.The implementation is based of 'Uniform Asymptotic Expansion of Bessel
    functions in terms of Airy functions' as described in 'Handbook of 
    Mathematical Functions' by Abramowitz and Stegun (1964).
  3.For the computation of complex valued Airy function, routines of Amos are 
    used. These are available at www.netlib.org  
  4.The emphasis is on the utility of the routines for large (complex) order
    and large argument. 
  5.Caution: Routines show numerical instabilities for the case when the 
    absolute value of ratio of order to argument is close to 1.  

  * Please report the bugs in the code to  k.r.hiremath@ieee.org

    Copyright (C) <2004>  <Kirankumar R. Hiremath>

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Library General Public
    License as published by the Free Software Foundation; either
    version 2 of the License, or (at your option) any later version.

    This library is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Library General Public License for more details.

    You should have received a copy of the GNU Library General Public
    License along with this library; if not, write to the
    Free Software Foundation, Inc., 59 Temple Place - Suite 330,
    Boston, MA  02111-1307, USA.

*/


/* Commnets:
   In BesselJ(), IERR = 11 means routine is not meant to handle situation 
   with order=0
   In DBesselY(), BesselY(), IERR = 12 means argument is zero, numerical 
   exception, division by zero.
 
  KODE = 1 => without scaling 
  KODE = 2 => with scaling    

  Error codes:
  IERR =  0 => :-) Success
  IERR = 11 => situation beyond scope of routines
  IERR = 12 => division by zero
  IERR =  5 => no computation, algorithm termination condition not met
  IERR =  4 => |z| too large. no computation.
  IERR =  3 => |z| too large. ans may wrong.
  IERR =  2 => overflow. no computation. Re(\zeta) too large for KODE = 1
  IERR =  1 => inout error

  
  Underflow indicator (in Airy function calculations): 
  NZ   = 0  => Normal return
  NZ   = 1   , AI=CMPLX(0.0D0,0.0D0) DUE TO UNDERFLOW IN
                              -PI/3.LT.ARG(Z).LT.PI/3 ON KODE=1
*/


#include "bessel.h"


const complex<double> I(0.0,1.0); // complex notation "i"

/* error message */
void complex_bessel_error(char *s)
{
	fprintf(stderr, " ERROR: complex_bessel: %s.\n", s);
	exit(1);
}
  
extern "C"
{
    extern void zairy_( double *zr, double *zi, int *ID, int *KODE, double *air, double *aii, int *NZ, int *IERR); //644 airy()
    extern void zbiry_( double *zr, double *zi, int *ID, int *KODE, double *air, double *aii, int *IERR); //644 biry()
    extern double dgamma_(double *x);
};

double lambda(int s);
double mu(int s);
complex<double> u(int k, complex<double> t);
complex<double> v(int k, complex<double> t);
complex<double> a(int k, complex<double> zeta, complex<double> z);
complex<double> b(int k, complex<double> zeta, complex<double> z);

complex<double> a(int k, complex<double> zeta); // new def
complex<double> b(int k, complex<double> zeta); // new def

complex<double> c(int k, complex<double> zeta, complex<double> z);
complex<double> d(int k, complex<double> zeta, complex<double> z);


//------------Function definitions start here-----------------


//-----------------arccos() starts-----------------------------
complex<double> arccos(complex<double> zin)
{  extern complex<double> I;
   return (log(zin + I*sqrt(1.0 - zin*zin))/I); 
 }
//------------------arccos() ends-----------------------------

 //-------------------Hankel function set start -------------------------
 complex<double> HankelH1(complex<double> order, complex<double> argument, int KODE, int *NZ, int *IERR)
{   
    extern complex<double> I;
    
    complex <double> J= BesselJ(order,argument, KODE, &(*NZ), &(*IERR));
    if((*IERR) > 10 )
    { 
	return complex<double> (nan, nan);
    }

    complex <double> Y= BesselY(order,argument, KODE, &(*NZ), &(*IERR));
    if((*IERR) > 10 )
    { 
	return complex<double> (nan, nan);
    }
    
     return ( J + I*Y ) ;
}

 complex<double> HankelH2(complex<double> order, complex<double> argument, int KODE, int *NZ, int *IERR)
 { 
     int NZ_tmp, IERR_tmp;

     extern complex<double> I;
     
     complex <double> J= BesselJ(order,argument, KODE, &NZ_tmp, &IERR_tmp);

     if((IERR_tmp) > 10 )
     { 
	 return complex<double> (nan, nan);
     }

     complex <double> Y= BesselY(order,argument, KODE, &NZ_tmp, &IERR_tmp);
     if((IERR_tmp) > 10 )
     { 
	 return complex<double> (nan, nan);
     }
    
     *NZ = NZ_tmp;
     *IERR = IERR_tmp;

     return ( J - I*Y ) ;
 }

 complex<double> DHankelH1(complex<double> order, complex<double> argument, int KODE, int *NZ, int *IERR)
 { 
     extern complex<double> I;
     
     complex <double> DJ= DBesselJ(order,argument, KODE, &(*NZ), &(*IERR));
     if((*IERR) > 10 )
     { 
	 return complex<double> (nan, nan);
     }

     complex <double> DY= DBesselY(order,argument, KODE, &(*NZ), &(*IERR));
     if((*IERR) > 10 )
     { 
	 return complex<double> (nan, nan);
     }
    
     return ( DJ + I*DY ) ;
 }


 complex<double> DHankelH2(complex<double> order, complex<double> argument, int KODE, int *NZ, int *IERR)
 { extern complex<double> I;

     complex <double> DJ= DBesselJ(order,argument, KODE, &(*NZ), &(*IERR));
     if((*IERR) > 10 )
     { 
	 return complex<double> (nan, nan);
     }

     complex <double> DY= DBesselY(order,argument, KODE, &(*NZ), &(*IERR));
     if((*IERR) > 10 )
     { 
	 return complex<double> (nan, nan);
     }

     return ( DJ - I*DY ) ;
 }

 //-----------------Hankel function set end-------------------------------


 //--------------DBesselJ() starts ------------------------------------
 complex<double> DBesselJ(complex<double> order, complex<double> argument, int KODE, int *NZ, int *IERR)
 {
     *NZ=0;
     *IERR=0;
     
   if (abs(argument) == 0)
     {
       return 0;
     }
   else if ( order == 0.0 )
    {
	*IERR = 11;
	return complex<double>(nan, nan);
    }
   else if (order == argument)
     {
       //using same notation of Abramowitz and Stegun
       double a = 0.4473073184;
       double gamma[4] = { 1.0, 0.0073015873, -0.0009373, 0.000444 };
       double b = 0.4108501939;
       double delta[4] = { 0.2, -947.0/346500.0, 0.0006047, -0.00038};

       complex<double> a1, a2, ans;


       a1  = gamma[0] + gamma[1]/pow(order,2.0) + gamma[2]/pow(order,4.0) 
	     +  gamma[3]/pow(order,6.0);
       a2 = delta[0] + delta[1]/pow(order,2.0) + delta[2]/pow(order,4.0) 
	     +  delta[3]/pow(order,6.0);
       ans = ((b / pow(order,2.0/3.0))*a1) - ((a/pow(order,4.0/3.0))*a2);

       return ans;
     }
   else
     {
       complex<double> z, zeta, ans;

       z = argument/order;

       if ( abs(z) < 1.0 )
       { 	zeta = pow(1.5*(log((1.0 + sqrt(1.0 - z*z))/z) -
				sqrt(1.0 - z*z)), 2.0/3.0);
       }
       else
       {
 	zeta = -1.0 * pow(1.5*(sqrt(z*z -1.0) - arccos(1.0/z)), 2.0/3.0);
       }

 //------------Computation of Ai and DAi------------------------------------

     complex<double> Airy_argument = pow(order, 2.0/3.0)*zeta;
     complex<double> AiryAi, DAiryAi, scale_factor;
     double zr, zi, air, aii;
     int Airy_ID;
     
     zr = real(Airy_argument);
     zi = imag(Airy_argument);

     //calling 644 routine to compute Ai(z)
     if (KODE == 1)
	{
	    Airy_ID=0; //calculating airy()
	    zairy_(&zr, &zi, &Airy_ID, &KODE, &air, &aii, NZ, IERR);
	    AiryAi=complex<double> (air, aii);

	    Airy_ID=1; //now calculating derivative of airy()
	    zairy_(&zr, &zi, &Airy_ID, &KODE, &air, &aii, NZ, IERR);
	    DAiryAi=complex<double>(air, aii);
	    
	}
	else if (KODE == 2)
	{
	    complex<double> scale_factor=exp(2.0*(Airy_argument*sqrt(Airy_argument))/3.0);
	    
	    Airy_ID=0; //calculating airy()
	    zairy_(&zr, &zi, &Airy_ID, &KODE, &air, &aii, NZ, IERR);
	    AiryAi=complex<double>(air, aii);
	    AiryAi =  AiryAi/scale_factor;

	    Airy_ID=1; //now calculating derivative of airy()
	    zairy_(&zr, &zi, &Airy_ID, &KODE, &air, &aii, NZ, IERR);
	    DAiryAi=complex<double>(air, aii);
	    DAiryAi = DAiryAi/scale_factor;
	    
	}
	else
	    cout << "Invalid KODE. Quitting. \n";

     //---------------------------------------------

     complex<double> a1, a2, a3, a4, a5;

     a1 = -2.0*sqrt(sqrt((1.0 - z*z)/(4.0*zeta)))/z;
     a2 = AiryAi/pow(order, 4.0/3.0);
     a3 = c(0,zeta,z); 
     a4 = DAiryAi/(pow(order, 2.0/3.0));
     a5 = d(0,zeta,z) + d(1, zeta,z)/(order*order);

     ans = a1*(a2*a3 + a4*a5);

     return ans  ;

     }
 }

 //--------------DBesselJ() finishs------------------------------------

 //--------------DBesselY() starts ------------------------------------
complex<double> DBesselY(complex<double> order, complex<double> argument, int KODE, int *NZ, int *IERR)
 {
     *NZ=0;
     *IERR=0;
          
 if (abs(argument) == 0)
     {
         //IERR or NZ has to be set appropiately :-)
	 *IERR=12;
         return complex<double>(nan,nan);
     }
    else if ( order == 0.0 )
    {
	*IERR = 11;
	return complex<double>(nan, nan);
    }
   else if (order == argument)
     {
       //using same notation of Abramowitz and Stegun
       double a = 0.4473073184;
       double gamma[4] = { 1.0, 0.0073015873, -0.000937300, 0.000444 };
       double b = 0.4108501939;
       double delta[4] = { 0.2, -947.0/346500.0, 0.0006047, -0.00038};

       complex<double> a1, a2, ans;


       a1  = gamma[0] + gamma[1]/pow(order,2) + gamma[2]/pow(order,4) +
 	    gamma[3]/pow(order,6);
       a2 = delta[0] + delta[1]/pow(order,2) + delta[2]/pow(order,4) +
 	          delta[3]/pow(order,6);
       ans = sqrt(3.0)*((b / pow(order,2.0/3.0))*a1) + ((a/pow(order,4.0/3.0))*a2);

       *IERR = 0;
       return ans;
     }
   else
     {
       complex<double> z, zeta, ans;
       z = argument/order;

       if ( abs(z) < 1.0 )
 	zeta = pow(1.5*(log((1.0 + sqrt(1.0 - z*z))/z) - sqrt(1.0 - z*z)), 2.0/3.0);
       else
 	zeta = -1.0 * pow(1.5*(sqrt(z*z -1.0) - arccos(1.0/z)), 2.0/3.0);


 //------------Computation of AiryBi and AiryDBi-------------------

     complex<double> Airy_argument = pow(order, 2.0/3.0)*zeta;
     complex<double> scale_factor, AiryBi, DAiryBi;
     double zr, zi, air, aii;
     int Airy_ID;
     
     zr = real(Airy_argument);
     zi = imag(Airy_argument);

     //calling 644 routine to compute Ai(z)
     if (KODE == 1)
	{
	    Airy_ID=0; //calculating airy()
	    zbiry_(&zr, &zi, &Airy_ID, &KODE, &air, &aii, IERR);
	    AiryBi=complex<double> (air, aii);

	    Airy_ID=1; //now calculating derivative of airy()
	    zbiry_(&zr, &zi, &Airy_ID, &KODE, &air, &aii, IERR);
	    DAiryBi=complex<double>(air, aii);
	    
	}
	else if (KODE == 2)
	{
	    complex<double> scale_factor=exp(fabs(real(2.0*(Airy_argument*sqrt(Airy_argument))/3.0)));
	    
	    Airy_ID=0; //calculating airy()
	    zbiry_(&zr, &zi, &Airy_ID, &KODE, &air, &aii, IERR);
	    AiryBi=complex<double>(air, aii);
	    AiryBi =  AiryBi*scale_factor;

	    Airy_ID=1; //now calculating derivative of airy()
	    zbiry_(&zr, &zi, &Airy_ID, &KODE, &air, &aii, IERR);
	    DAiryBi=complex<double>(air, aii);
	    DAiryBi = DAiryBi*scale_factor;

	}
	else
	    cout << "Invalid KODE. Quitting. \n";
	


 //---------------------------------------------

     complex<double> a1, a2, a3, a4, a5;

     a1 = 2.0*pow((1.0 - z*z)/(4.0*zeta), 0.25)/z;
     a2 = AiryBi/pow(order, 4.0/3.0);
     a3 = c(0,zeta,z); 
     a4 = DAiryBi/pow(order, 2.0/3.0);
     a5 = d(0,zeta,z) + d(1, zeta,z)/pow(order,2);

     ans =a1*(a2*a3 + a4*a5);

     return ans  ;

     }
 }


 //--------------DBesselY() finishs------------------------------------


//--------------------------BesselY() starts--------------------------

complex<double> BesselY(complex<double> order, complex<double> argument, int KODE, int *NZ, int *IERR)
{
    //still NZ and IERR has to be used systematically :-)
    *NZ=0;
    *IERR=0; *IERR = 0;

    if (abs(argument) == 0)
    {
	*IERR=12;
	return complex<double>(nan,nan);
    }
    else if ( order == 0.0 )
    {
	*IERR = 11;
	return complex<double>(nan, nan);
    }
    else if (order == argument)
    {
	//using same notation of Abramowitz and Stegun
	double a = 0.4473073184f;
	double alpha[4] = { 1.0, -1.0/225.0, 0.000693735, -0.0003538 };
	double b = 0.4108501939f;
	double beta[4] = {1.0/70.0, -1213.0/1023750.0, 0.0004378, -0.00038};
	
	complex<double> a1, a2, ans;


	a1  = alpha[0] + alpha[1]/pow(order,2) + alpha[2]/pow(order,4) + 
              alpha[3]/pow(order,6);
	a2 = beta[0] + beta[1]/pow(order,2) + beta[2]/pow(order,4) + 
             beta[3]/pow(order,6);
	ans = -(((sqrt(3.0)*a / pow(order,1.0/3.0))*a1) + ((sqrt(3.0)*b/pow(order,5.0/3.0))*a2));


	 *IERR = 0;
	return ans;
    }
    else
    {
	complex<double> z, zeta, ans;

	z = argument/order;

	if ( abs(z) < 1.0 )
	    zeta = pow(1.5*(log((1.0 + sqrt(1.0 - z*z))/z) - sqrt(1.0 - z*z)), 2.0/3.0);
	else
	    zeta = -1.0 * pow(1.5*(sqrt(z*z -1.0) - arccos(1.0/z)), 2.0/3.0);


	//------------Computation of AiryBi and DAiryBi--------------------

	complex<double> Airy_argument = pow(order, 2.0/3.0)*zeta;
	complex<double> AiryBi, DAiryBi, scale_factor;
	double zr, zi, air, aii;
	int Airy_ID;
	zr=real(Airy_argument);
	zi=imag(Airy_argument);
	
	
    
	if (KODE == 1)
	{
	    Airy_ID=0; //calculating biry()
	    zbiry_(&zr, &zi, &Airy_ID, &KODE, &air, &aii, IERR);
	    AiryBi=complex<double> (air, aii);

	    Airy_ID=1; //now calculating derivative of airy()
	    zbiry_(&zr, &zi, &Airy_ID, &KODE, &air, &aii, IERR);
	    DAiryBi=complex<double>(air, aii);
	    
	}
	else if (KODE == 2)
	{
	    complex<double> scale_factor=exp(fabs(real(2.0*(Airy_argument*sqrt(Airy_argument))/3.0)));
	    
	    Airy_ID=0; //calculating airy()
	    zbiry_(&zr, &zi, &Airy_ID, &KODE, &air, &aii, IERR);
	    AiryBi=complex<double>(air, aii);
	    AiryBi =  AiryBi*scale_factor;

	    Airy_ID=1; //now calculating derivative of airy()
	    zbiry_(&zr, &zi, &Airy_ID, &KODE, &air, &aii, IERR);
	    DAiryBi=complex<double>(air, aii);
	    DAiryBi = DAiryBi*scale_factor;
	    

	}
	else
	    cout << "Invalid KODE. Quitting. \n";
	
	

	complex<double> a1, a2, a3, a4, a5;
	a1 = sqrt(sqrt(4.0*zeta/(1.0 - z*z)));
	a2 = AiryBi/pow(order, 1.0/3.0);
	a3 = a(0,zeta) +a(1, zeta)/(order*order); 
	a4 = DAiryBi/pow(order, 5.0/3.0);
	a5 = b(0,zeta) + b(1, zeta)/(order*order);  

	ans =-a1*(a2*a3 + a4*a5);

	return ans  ;
	
    }
}


//--------------------------BesselY() ends--------------------------  
 
//--------------------------BesselJ() starts-------------------------- 

complex<double> BesselJ(complex<double> order, complex<double> argument, int KODE, int *NZ, int *IERR)
{
    //still NZ and IERR has to be used systematically :-)
    *NZ=0;
    *IERR=0;

    if (abs(argument) == 0)
    {
 	return 0;
    }
    else if ( order == 0.0 )
    {
	*IERR = 11;
	return complex<double>(nan, nan);
    }
    else if (order == argument)
    {
	//using same notation of Abramowitz and Stegun
	double a = 0.4473073184;
	double alpha[4] = { 1.0, -1.0/225.0, 0.000693735, -0.0003538 };
	double b = 0.4108501939;
	double beta[4] = {1.0/70.0, -1213.0/1023750.0, 0.0004378, -0.00038};
	
	complex<double> a1, a2, ans;


	a1  = alpha[0] + alpha[1]/pow(order,2) + alpha[2]/pow(order,4) +  
              alpha[3]/pow(order,6);
	a2 = beta[0] + beta[1]/pow(order,2) + beta[2]/pow(order,4) +  
             beta[3]/pow(order,6);
	ans = ((a / pow(order,1.0/3.0))*a1) - ((b/pow(order,5.0/3.0))*a2);

	 *IERR = 0;
	return ans;
    }
    else
    {
	complex<double> z, zeta, ans;

	z = argument/order;

	if ( abs(z) < 1.0)
	    zeta = pow(1.5*(log((1.0 + sqrt(1.0 - z*z))/z) - sqrt(1.0 - z*z)), 2.0/3.0);
	else
	    zeta = -1.0 * pow(1.5*(sqrt(z*z -1.0) - arccos(1.0/z)), 2.0/3.0);


	//------------Computation of AiryAi and DAiryAi--------------------

	complex<double> Airy_argument = pow(order, 2.0/3.0)*zeta;
	
	complex<double> AiryAi, DAiryAi, scale_factor;
	double zr, zi, air, aii;
	int Airy_ID;
	zr=real(Airy_argument);
	zi=imag(Airy_argument);
	
	
    
	if (KODE == 1)
	{
		Airy_ID=0; //calculating airy()
		zairy_(&zr, &zi, &Airy_ID, &KODE, &air, &aii, NZ, IERR);
		AiryAi=complex<double> (air, aii);

		Airy_ID=1; //now calculating derivative of airy()
		zairy_(&zr, &zi, &Airy_ID, &KODE, &air, &aii, NZ, IERR);
		DAiryAi=complex<double>(air, aii);

	}
	else if (KODE == 2)
	{
	    scale_factor=exp(2.0*(Airy_argument*sqrt(Airy_argument))/3.0);

	    Airy_ID=0; //calculating airy()
	    zairy_(&zr, &zi, &Airy_ID, &KODE, &air, &aii, NZ, IERR);
	    AiryAi=complex<double>(air, aii);
	    AiryAi =  AiryAi/scale_factor;

	    Airy_ID=1; //now calculating derivative of airy()
	    zairy_(&zr, &zi, &Airy_ID, &KODE, &air, &aii, NZ, IERR);
	    DAiryAi=complex<double>(air, aii);
	    DAiryAi = DAiryAi/scale_factor;
	    
	}
	else
	    cout << "Invalid KODE. Quitting. \n";
	
	

	complex<double> a1, a2, a3, a4, a5;
	a1 = sqrt(sqrt(4.0*zeta/(1.0 - z*z)));
	a2 = AiryAi/pow(order, 1.0/3.0);
	a3 = a(0,zeta) + a(1,zeta)/(order*order); 
	a4 = DAiryAi/pow(order, 5.0/3.0);
	a5 = b(0,zeta)  + b(1, zeta)/(order*order); 

	ans =a1*(a2*a3 + a4*a5);

	return ans  ;
	
    }
}

//---------------------------BesselJ() finish---------------


//----------------lambda()------------------------------
double lambda(int s)
{
    int i, ifail;
    double out_lambda=1.0, s_fact, s_in_double;
    s_in_double = s + 1.0;
    
    
    for (i=0; i<= (2*s -1); i++)
    {
 	out_lambda = out_lambda*(2*s + (2*i+1));
     }
    
    s_fact=dgamma_(&s_in_double);

    out_lambda = out_lambda /(s_fact * pow(144.0,s));
    return out_lambda;
    
}

//----------------mu()---------------------------------
double mu(int s)
{ return -1.0*(6.0*s + 1.0)*lambda(s)/(6.0*s -1.0);}


//---------------u()-----------------------------------

complex<double> u(int k, complex<double> t)
{
    if (k == 0 )  return 1.0;
    else if ( k == 1) return (3.0*t - 5.0*pow(t,3.0))/24.0;
    else if ( k == 2)return  (81.0*pow(t,2.0) - 462.0*pow(t,4.0) + 385.0*pow(t,6.0))/1152.0;
    else if ( k == 3) return (30375.0*pow(t,3) - 369603.0*pow(t,5.0) + 765765.0*pow(t,7.0) - 425425.0*pow(t,9.0))/414720.0;
    else  if ( k == 4) return (4465125.0*pow(t,4.0) - 94121676.0*pow(t,6.0) + 349922430.0*pow(t,8.0) - 446185740.0*pow(t,10.0) + 185910725.0*pow(t,12.0))/39813120.0;
    else if ( k == 5) return -pow(t,5.0)/6688604160.0*(-1519035525.0 + 49286948607.0*pow(t,2.0) - 284499769554.0*pow(t,4.0) + 614135872350.0*pow(t,6.0) - 566098157625.0*pow(t,8.0) + 188699385875.0*pow(t,10.0));
    else if (k == 6) return pow(t,6.0)/4815794995200.0*(2757049477875.0 - 127577298354750.0 * pow(t,2.0) + 1050760774457901.0*pow(t,4.0) - 3369032068261860.0*pow(t,6.0) + 5104696716244125.0*pow(t,8.0) - 3685299006138750.0*pow(t,10.0) + 1023694168371875.0*pow(t,12.0));
    //consider out ouf range case
    
}

//---------------v()-----------------------------------

complex<double> v(int k, complex<double> t)
{
    if (k == 0 )  return 1;
    else if ( k == 1) return (-9.0*t + 7.0*pow(t,3.0))/24.0;
    else if ( k == 2)return  (-135.0*pow(t,2.0) + 594.0*pow(t,4.0) - 455.0*pow(t,6.0))/1152.0;
    else if ( k == 3) return (-42525.0*pow(t,3.0) + 451737.0*pow(t,5.0) -  883575.0*pow(t,7.0) + 475475.0*pow(t,9.0))/414720.0;
    
    //consider out of range case

}

//-------------a()--------------------------------------
complex<double> a(int k, complex<double> zeta)
{
    complex<double> eta = zeta / pow(2.0, 1.0/3.0);
    
    if ( k == 0)
	return complex<double>(1.0, 0);
    else if (k ==1 )
	return -1.0/225.0 - 71.0 * eta / 38500.0 + 82.0*eta*eta/73125.0;
    else if (k == 2)
	return 151439.0/218295000.0 + 68401.0 * eta / 147262500.0 - 1796498167.0*eta*eta/4193689500000.0 - 583721053.0*eta*eta*eta/830718281250.0;
    else if (k == 3)
	return - 887278009.0/2504935125000.0 - 3032321618951.0*eta/9708942993750000.0;
    else
	{
	    cout <<"Error : Code is not implemented for k > 3. \n";
	    return complex<double> (nan, nan);
	}
}

//-------------b()--------------------------------------
complex<double> b(int k, complex<double> zeta)
{
    complex<double> eta = zeta / pow(2.0, 1.0/3.0);
    
    if ( k == 0)
	return pow(2.0, 1.0/3.0)*(1.0/70.0 + 2.0 *eta/225.0 + 138.0*eta*eta/67373.0 - 296.0*eta*eta*eta/511875.0 - 38464.0*eta*eta*eta*eta/63669375.0);
    else if (k ==1 )
	return pow(2.0, 1.0/3.0)*(-1213.0/1023750.0 -3757.0*eta/2695000.0 - 3225661.0*eta*eta/6700443750.0 + 90454643.0*eta*eta*eta/336992906250.0);
    else if (k == 2)
	return pow(2.0, 1.0/3.0)*(16542537833.0/37743205500000.0 + 115773498223.0*eta/162820783125000.0 + 548511920915149.0*eta*eta/1721719224225000000.0);
    else if (k == 3)
	return pow(2.0, 1.0/3.0)*(-430990563936859253.0/568167343994250000000.0 - 3191320338955050557.0*eta/7777535495585625000000.0);
    else
	{
	    cout <<"Error : Code is not implemented for k > 3. \n";
	    return complex<double> (nan, nan);
	}
}

//-------------c()--------------------------------------
complex<double> c(int k, complex<double> zeta, complex<double> z)
{
    
    int s;
    complex<double> out_a=(0,0);
    for(s=0; s<= 2*k+1; s++)
    {
	out_a = out_a + pow(zeta, -1.5*s) * v(2*k -s +1, pow(1.0-z*z, -0.5))*mu(s);
    }
    
    return -1.0*sqrt(zeta)*out_a;
}

//-------------d()--------------------------------------
complex<double> d(int k, complex<double> zeta, complex<double> z)
{
    int s;
    complex<double> out_b=(0,0);
    for(s=0; s<= 2*k; s++)
    {   double lam = lambda(s);

    out_b = out_b + lam*(pow(zeta, -1.5*s) * v(2*k -s, 1.0/sqrt(1.0-z*z)));
    
    }
    
    return out_b;
}
